vec4 cubicFilter(sampler3D volume, in VOLUME_STRUCT volumeStruct, in vec3 sampleOrg, float B, float C) {
	const int radius = 2;
	vec4 result = getVoxel(volume, volumeStruct, sampleOrg);
	if (radius < 1) return result;
	
	vec3 samplePoint = sampleOrg * volumeStruct.datasetDimensions_;
	vec3 pivotPoint = round(samplePoint);
	vec3 distance = samplePoint - pivotPoint;
	const float eps = 0.001;
	if (length(distance) < eps) return result;

	result.a = 0.0;
	const int width = 2 * radius + 1;
	
	int pos;
	vec3 arg;
	vec3 argSqr;
	vec3 argCube;
	vec3 spline[width];
	for (int i = -radius; i <= radius; i++) {
		pos = i + radius;
		arg = abs(distance - i);
		argSqr = arg * arg;
		argCube = arg * argSqr; 
		for (int dim = 0; dim < 3; dim++) {
			if (arg[dim] < 1) {
				spline[pos][dim] = ((12 - 9*B - 6*C) * argCube[dim] + (-18 + 12*B + 6*C) * argSqr[dim] + (6 - 2*B)) / 6;
			} 
			else {
				if (arg[dim] < 2) {
					spline[pos][dim] = ((-B - 6*C) * argCube[dim] + (6*B + 30*C) * argSqr[dim] + (-12*B - 48*C) * arg[dim] + (8*B + 24*C)) / 6;
				}
				else {
					spline[pos][dim] = 0.0;
				}
			}
		}
	}

	// convolution-sum (3D)
	for (int x = -radius; x <= radius; x++) {
        for (int y = -radius; y <= radius; y++) {
            for (int z = -radius; z <= radius; z++) {
				result.a += getVoxel(volume, volumeStruct, (pivotPoint + vec3(x, y, z)) / volumeStruct.datasetDimensions_).a  
							* spline[x + radius].x * spline[y + radius].y * spline[z + radius].z;
			}
		}
	}

	return result;
}